Podsumowanie

(poniższe komentarze dostępne są również w miejscu ich dotyczącym)

a.d.3. Kod pozwalający wczytać dane z pliku.

Do wczytania danych wykorzystałem funkcję fread z pakietu data.table. Funkcja ta jest znacznie szybsza niż read.table i pozwala na usuwanie zbędnych kolumn już na etapie wczytywania danych. Dodatkowo prawidłowo rozpoznaje kolumny o typie Integer

a.d.5. Kod przetwarzający brakujące dane.

Dane zostają podzielone na 3 części (part_00, part_01 i part_02). Prefisky kolumn usunięte, dodana kolumna part z wartością odpowiedniego partu będącego źródłem danych. Następnie dane są łączone w jeden data frame. Proces ten usuwa wszystkie wartości NA - uzyskane dane są więc wyczyszczone.

Jeżeli any_meas* = 0, to zakładam, że nie nastąpił‚ pomiar (w kolejnych kolumnach dane NA, gdyż nie było podstawy do obliczeń). W takiej sytuacji zamiana NA na np. średnią nie ma sensu (0 elektronów nie może mieć kształtu itp.) dlatego też dane dotyczące danego part zostają usunięte.

*any_meas = shape_segments_count + density_segments_count + volume + electrons + mean + std + max + max_over_std + skewness + parts) [żadna z tych zmiennych nie przyjmuje wartości < 0, więc suma == 0 jednoznacznie wskazuje na brak tych danych]

a.d.8. Sekcję sprawdzającą korelacje między zmiennymi

Poniżej znajduje się wykres korelacji między zmiennymi. Zauważyć można, że:

Duża część zmiennych jest w silnej korelacji z innymi zmiennymi
Zmienne density_Z* i shape_z* są w bardzo silnej korelacji (duże czerwone prostokąty)
Zmienne local* korelują z zmiennymi z part (były obliczane na ich podstawie)

a.d.14 i 15

Spośród testowanych metod najlepiej sprawdziły się:

Linear Regression (lm) dla regresji
Random Forest (rf) dla klasyfikacji

Z racji dużego zbioru danych (ok. 1 mln rekordów) za metodę wybierania zbioru walidacyjnego wybrałem bootstrap (losowanie ze zwracaniem) - metoda zapewnia stratyfikację danych (https://machinelearningmastery.com/how-to-estimate-model-accuracy-in-r-using-the-caret-package/)

a.d.14.Sekcję sprawdzającą czy na podstawie wartości innych kolumn można przewidzieć:

Na podstawie obliczeń usatliłem, że odpowiednio 3% i 1% danych powinien wystarczyć jako zbiór testowy (z obliczeniami można zapoznać się w sekcji 14)

Liczba elektronów została określona z miarami:

RMSE = 64,453

r^2 = 0,48

Liczba elektronów została określona z miarami:

RMSE = 9,45

r^2 = 0,49

a.d.15.Sekcję próbującą stworzyć klasyfikator przewidujący wartość atrybutu res_name

Na podstawie obliczeń usatliłem, że 1% danych powinien wystarczyć jako zbiór testowy (z obliczeniami można zapoznać się w sekcji 15)

Random Forest uzyskał precyzję = 0,6

Klasyfikator naiwny (wskazujący najliczniejszą klasę) dokonałby predykcji res_name = SO4. Miara Accuracy = 0,15.

1. Kod wyliczający wykorzystane biblioteki.

library(data.table)
library(DT)
library(ggplot2)
library(plotly)
library(dplyr)
library(reshape)
library(caret)

2. Kod zapewniający powtarzalność wyników.

set.seed(123)

3. Kod pozwalający wczytać dane z pliku.

out_columns <- c("blob_coverage", "res_coverage", "title", "pdb_code", "res_id", "chain_id", "blob_volume_coverage", "blob_volume_coverage_second", "res_volume_coverage", "res_volume_coverage_second", "skeleton_cycle_4", "skeleton_diameter", "skeleton_cycle_6", "skeleton_cycle_7", "skeleton_closeness_006_008", "skeleton_closeness_002_004", "skeleton_cycle_3", "skeleton_avg_degree", "skeleton_closeness_004_006", "skeleton_closeness_010_012", "skeleton_closeness_012_014", "skeleton_edges", "skeleton_radius", "skeleton_cycle_8_plus", "skeleton_closeness_020_030", "skeleton_deg_5_plus", "skeleton_closeness_016_018", "skeleton_closeness_008_010", "skeleton_closeness_018_020", "skeleton_average_clustering", "skeleton_closeness_040_050", "skeleton_closeness_014_016", "skeleton_center", "skeleton_closeness_000_002", "skeleton_density", "skeleton_closeness_030_040", "skeleton_deg_4", "skeleton_deg_0", "skeleton_deg_1", "skeleton_deg_2", "skeleton_deg_3", "skeleton_graph_clique_number", "skeleton_nodes", "skeleton_cycles", "skeleton_cycle_5", "skeleton_closeness_050_plus", "skeleton_periphery", "local_cut_by_mainchain_volume", "local_near_cut_count_C", "local_near_cut_count_other", "local_near_cut_count_S", "local_near_cut_count_O", "local_near_cut_count_N", "fo_col", "fc_col", "weight_col", "grid_space", "solvent_radius", "solvent_opening_radius", "resolution_max_limit", "part_step_FoFc_std_min", "part_step_FoFc_std_max", "part_step_FoFc_std_step", "skeleton_data", "local_res_atom_count", "local_res_atom_non_h_occupancy_sum", "local_res_atom_non_h_electron_occupancy_sum", "local_res_atom_C_count", "local_res_atom_N_count", "local_res_atom_O_count", "local_res_atom_S_count", "dict_atom_C_count", "dict_atom_N_count", "dict_atom_O_count", "dict_atom_S_count")

file = "C:/Users/Developer/Desktop/all_summary.csv"
data <- fread(file = file, 
             sep = ";", 
             header = TRUE,
             na.string = c(",,", "NAN", "nan"),
             drop = out_columns
            )

4. Kod usuwający z danych wiersze posiadające wartość zmiennej res_name

out_res_names = c("UNK", "UNX", "UNL", "DUM", "N", "BLOB", "ALA", "ARG", "ASN", "ASP", "CYS", "GLN", "GLU", "GLY", "HIS", "ILE", "LEU", "LYS", "MET", "MSE", "PHE", "PRO", "SEC", "SER", "THR", "TRP", "TYR", "VAL", "DA", "DG", "DT", "DC", "DU", "A", "G", "T", "C", "U", "HOH", "H20", "WAT")
data <- data %>% filter(!res_name %in% out_res_names)

5. Kod przetwarzający brakujące dane.

data_part00 <- data %>% 
    select(-starts_with("part_01"), -starts_with("part_02")) %>% 
    rename_at(.vars = vars(starts_with("part_00_")), .funs = funs(sub("^part_00_", "", .))) %>%
    mutate(any_meas = shape_segments_count + density_segments_count + volume + electrons + mean + std + max + max_over_std + skewness + parts) %>%       
    filter(any_meas != 0) %>% 
    select(-any_meas) %>%
    mutate(part = "part_00")

data_part01 <- data %>% 
    select(-starts_with("part_00"), -starts_with("part_02")) %>% 
    rename_at(.vars = vars(starts_with("part_01_")), .funs = funs(sub("^part_01_", "", .))) %>%
    mutate(any_meas = shape_segments_count + density_segments_count + volume + electrons + mean + std + max + max_over_std + skewness + parts) %>%       
    filter(any_meas != 0) %>% 
    select(-any_meas) %>%  
    mutate(part = "part_01")

data_part02 <- data %>% 
    select(-starts_with("part_00"), -starts_with("part_01")) %>% 
    rename_at(.vars = vars(starts_with("part_02_")), .funs = funs(sub("^part_02_", "", .))) %>%
    mutate(any_meas = shape_segments_count + density_segments_count + volume + electrons + mean + std + max + max_over_std + skewness + parts) %>%       
    filter(any_meas != 0) %>% 
    select(-any_meas) %>%  
    mutate(part = "part_02")
data <- rbind(data_part00, data_part01, data_part02)

W kolumnach znajduje się local_min - kolumna ta ma 0 we wszystkich wierszach - jest więc zbędna

data %>% select(local_min) %>% distinct()

##   local_min
## 1         0

data <- data %>% select(-local_min)
rm(data_part00, data_part01, data_part02)

6. Sekcję podsumowującą rozmiar zbioru i podstawowe statystyki.

Wymiary:

dim(data)

## [1] 1710450     125

Liczba unikalnych res_name

data %>% select(res_name) %>% distinct() %>% nrow()

## [1] 19603

Przykładów/part

data %>% group_by(part) %>% summarize(n = n()) %>% arrange(desc(n)) %>% prettyTable()

Podstawowe statystyki pozostałych kolumn

data %>% select(-res_name, -starts_with("part")) %>%
  summary() %>% 
  unclass() %>%
  data.frame(check.names = FALSE, stringsAsFactors = FALSE) %>% 
  prettyTable()

7. Kod ograniczający liczbę klas (res_name) do 50 najpopularniejszych wartości.

res_name50 <- data %>% group_by(res_name) %>% summarize(n = n()) %>% arrange(desc(n)) %>% head(50)
data <- data %>% filter(res_name %in% res_name50$res_name)
rm(res_name50)

8. Sekcję sprawdzającą korelacje między zmiennymi

Poniżej znajduje się wykres korelacji między zmiennymi. Zauważyć można, że: * Duża część zmiennych jest w silnej korelacji z innymi zmiennymi * Zmienne density_Z* i shape_z* są w bardzo silnej korelacji (duże czerwone prostokąty) * Zmienne local* korelują z zmiennymi z part (były obliczane na ich podstawie)

melted <- data %>% select(-res_name, -part) %>% cor() %>% melt
breaks <- sort(colnames(data))[seq(1, ncol(data), by = 6)]

(ggplot(data = melted, aes(x=X1, y=X2, fill=value)) + 
    theme(
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(angle=45)
      ) + 
    scale_x_discrete(breaks = breaks) + 
    scale_y_discrete(breaks = breaks) + 
    geom_tile() +
    scale_fill_gradient2(
        low = "blue", high = "red", mid = "white", 
        midpoint = 0, limit = c(-1,1))
      ) %>%
   ggplotly()

rm(melted, breaks)

9. Określenie ile przykładów ma każda z klas.

data %>% select(res_name) %>% group_by(res_name) %>% summarize(n = n()) %>% arrange(desc(n)) %>% prettyTable()

10. Wykresy rozkładów liczby atomów i elektronów.

(ggplot(data, aes(local_res_atom_non_h_count, fill = "red")) + geom_bar() + theme_bw()) %>% ggplotly()

(ggplot(data, aes(local_res_atom_non_h_electron_sum, fill = "red")) + geom_bar(width = 3) + theme_bw()) %>% ggplotly()

11. Tabelę pokazującą 10 klas z największą niezgodnością liczby atomów i liczby elektronów

diff <- data %>% select(res_name, local_res_atom_non_h_count, dict_atom_non_h_count, local_res_atom_non_h_electron_sum, dict_atom_non_h_electron_sum) %>%
  group_by(res_name) %>% summarize(atom_diff = sum(abs(local_res_atom_non_h_count - dict_atom_non_h_count)), electron_diff = sum(abs(local_res_atom_non_h_electron_sum - dict_atom_non_h_electron_sum)))

diff %>% select(-electron_diff) %>% arrange(desc(atom_diff)) %>% head(10) %>% prettyTable()

diff %>% select(-atom_diff) %>% arrange(desc(electron_diff)) %>% head(10) %>% prettyTable()

rm(diff)

12. Sekcję pokazującą rozkład wartości part_01

melted <- data %>% filter(part == "part_01") %>% select(shape_segments_count:density_Z_4_0) %>% melt
means <-  melted %>% group_by(variable) %>% summarise(mean=mean(value)) 
melted %>% ggplot(aes(value)) +
    geom_density() +
    geom_vline(data = means, aes(xintercept=mean), linetype="dashed", color = "red") +
    geom_text(data=means, mapping=aes(x=mean, y=0, label=signif(mean, digits = 4)), 
              size=3, angle=90, vjust= 1, hjust=0, color = "red"
      )+  
    facet_wrap(~variable, ncol = 4, scales = "free") +
    theme_bw()+
    theme(
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(angle=90)
      )

rm(melted, means)

13. Interaktywny wykres lub animację.

W pkt 8. i 10.

14.Sekcję sprawdzającą czy na podstawie wartości innych kolumn można przewidzieć:

liczbę elektronów

Posiadam ok 1mln przykładów podzielonych na n klas

n = data %>% select(local_res_atom_non_h_electron_sum) %>% distinct() %>% count()
n

## # A tibble: 1 x 1
##       n
##   <int>
## 1   176

Zbiór testowy złożony z 3% tych danych powinien być wystarczający. Każda z klas będzie miała więc ok. k przykładów (zależnie od rozkładu)

k = nrow(data) * 0.03 / n
k

##          n
## 1 190.3739

idx <- createDataPartition(data$local_res_atom_non_h_electron_sum, p=0.03, list=F, times = 1)
test <- data.frame(data[idx,]) %>% select(-res_name, -local_res_atom_non_h_count)
train <- data.frame(data[-idx,]) %>% select(-res_name, -local_res_atom_non_h_count)

ctrl <- trainControl(
    method = "boot", number = 20)
 
fit <- train(local_res_atom_non_h_electron_sum ~ .,
             data = train,
             method = "lm",
             trControl = ctrl,
             metric = "Rsquared",
             maximize = TRUE)

fit

## Linear Regression 
## 
## 1083352 samples
##     120 predictor
## 
## No pre-processing
## Resampling: Bootstrapped (20 reps) 
## Summary of sample sizes: 1083352, 1083352, 1083352, 1083352, 1083352, 1083352, ... 
## Resampling results:
## 
##   RMSE      Rsquared   MAE     
##   64.43256  0.4826822  36.23728
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

rfClasses <- predict(fit, newdata = test) %>% round %>% as.integer

rm(idx, train, test)

liczbę atomów